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ABSTRACT 

We describe a Novel form of Adaptive softening (NovA) for collisionless A-body 
simulations, implemented in the RAMSES adaptive mesh refinement code. In RAM¬ 
SES - that we refer to as a ‘standard A-body method’ - cells are only split if they 
contain more than eight particles (a mass refinement criterion). Here, we introduce 
an additional criterion that the particle distribution within each cell be sufficiently 
isotropic, as measured by the ratio of the maximum to minimum eigenvalues of its 
moment of inertia tensor: rj = Amax/Amin- In this way, collapse is only refined if it 
occurs along all three axes, ensuring that the softening e is always of order twice the 
largest inter-particle spacing in a cell. This more conservative force softening crite¬ 
rion is designed to minimise spurious two-body effects, while maintaining high force 
resolution in collapsed regions of the flow. 

We test Nov A using an antisymmetric perturbed plane wave collapse (‘Valinia’ 
test) before applying it to warm dark matter (WDM) simulations. For the Valinia test, 
we show that - unlike the standard A-body method - NovA produces no numerical 
fragmentation while still being able to correctly capture fine caustics and shells around 
the collapsing regions. For the WDM simulations, we find that NovA converges sig¬ 
nificantly more rapidly than standard A-body, producing little or no spurious halos on 
small scales. We show, however, that determining whether or not halos exist below the 
free streaming mass Mfg is complicated by the fact that our halo finder (AHF) likely 
incorrectly labels some caustics and criss-crossing filaments as halos, while one or two 
particularly massive filaments appear to fragment in any version of NovA where re¬ 
finement is allowed. Such massive filaments may be physically unstable to collapse, as 
is the case for infinite, static, self-gravitating cylinders. We will use NovA in forth¬ 
coming papers to study the issue of halo formation below Mfg; filament stability; and 
to obtain new constraints on the temperature of dark matter. 

Key words: 


1 INTRODUCTION 

The A-body method is widely used for modelling the non¬ 
linear growth of structure in the Universe (e.g. Dehnen 
V Read, 2011; Kuhlen et ah, 2012). For collisionless non- 
relativistic (‘cold’) dark matter, it has been shown to be re¬ 
markably accurate, producing robust results that are numer¬ 
ically well-converged across different implementations (e.g. 
Heitmann et ah, 2008; Stadel et ah, 2009; Springel et ah, 
2008; Kim et ah, 2014). Such simulations provide an excel- 
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lent match to the observed large scale structure in the Uni¬ 
verse (e.g. Springel et ah, 2006), though on smaller scales - 
where baryons likely play a role (e.g. Navarro et ah, 1996; 
Read V Gilmore, 2005; Mashchenko et ah, 2008; Pontzen 
V Governato, 2011) - there are known discrepancies (e.g. 
Flores V Primack, 1994; Moore, 1994; Moore et ah, 1999; 
Klypin et ah, 1999). 

Despite the successes of the A-body method, since its 
inception there have been concerns about the effect of dis¬ 
creteness errors on numerical accuracy and convergence (e.g. 
Splinter et ah, 1998; Melott et ah, 1997; Diemand et ah, 
2004; Binney, 2004; Wang V White, 2007; Romeo et ah. 
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2008; Joyce et aL, 2009). These arise because the dark mat¬ 
ter fluid is represented by a set of ‘particles’, each with mass 
typically in the range 10^ — 10® Mq. To avoid spurious scat¬ 
tering between these particles the force is softened, for ex¬ 
ample using ‘Plummer’ (Plummer, 1915) softening: 

(e2 + |x« - x,f )3/2 

where Yij is the force between two particles i and j at posi¬ 
tions G is Newton’s gravitational constant; and e is the 
force softening. Equation 1 ensures that the force is clipped 
at a constant value as two particles approach one another. 
This reduces spurious two-body scattering, but it does not 
prevent numerical relaxation from occurring; that can only 
be combated by raising the number of particles (e.g. Dehnen, 
2001; Power et ah, 2003a; Binney, 2004; Diemand et ah, 
2004; Dehnen Sz Read, 2011). 

For cold dark matter (CDM) simulations, numerical 
convergence appears to be very good (e.g. Heitmann et ah, 
2008). However, discreteness errors may yet play a role when 
attempting to calculate power spectra, mass functions and 
higher order halo statistics at percent level accuracy, as will 
be required by next generation cosmological probes (e.g. 
Reed et ah, 2013; Amendola et ah, 2013). More problem¬ 
atic, however, are simulations that model a sharp cut-off in 
the initial power spectrum, as in warm dark matter^ (WDM; 
Bode et al. 2001a; Avila-Reese et al. 2001), or exotic infla¬ 
tionary models (Zentner & Bullock, 2003). The first WDM 
simulations appeared to find evidence of fragmentation - 
smaller halos forming later than larger ones - as evidenced 
by a sharp upturn in the halo mass function (e.g. Bode et ah, 
2001a; Avila-Reese et ah, 2001). However, this has now been 
traced to the numerical fragmentation of filaments due to 
discreteness effects. This is particularly worrisome since the 
‘spurious halos’ that form via this process diminish with par¬ 
ticle number only as leading to extremely slow con¬ 

vergence (Wang & White, 2007). 

The likely reason for the formation of ‘spurious halos’ 
in WDM simulations was only recently elucidated by Hahn 
et al. (2013). Using a new method for evolving collisionless 
fluids - where they track the foliations of the the 3D dark 
matter phase sheet - they find that the spurious halos result 

^ In WDM, it is supposed that the dark matter is non-relativistic 
for a time after decoupling, leading to a suppression in the growth 
of structure on small scales and at early times (e.g. Bode et al., 
2001a; Avila-Reese et al., 2001). Typically, this is modelled as 
an exponential cut-off in the initial power spectrum and indeed 
throughout this paper, where we refer to WDM simulations, this 
is what we mean. Fully self-consistent WDM models (for exam¬ 
ple, sterile neutrinos) have more complex model-dependent power 
spectra than this (e.g. Boyarsky et al., 2009). Furthermore, for 
hot dark matter, it can also become important to model the pri¬ 
mordial velocity dispersion of the dark matter particles. This has 
been attempted only a few times in the literature, most likely 
because of the computational cost involved. A proper treatment 
requires us to replace each ‘cold dark matter’ particle in the ini¬ 
tial conditions by ~ 1000 — 10, 000 lighter particles in order to 
well-sample the local velocity distribution function at each point 
in the flow (e.g. Avila-Reese et al., 2001; Hahn et al., 2013). As far 
as the authors are aware, at the time of writing, such an expensive 
approach has never been attempted. 


from large anisotropic force errors. With a more accurate 
force (as calculated by their new method), the spurious halos 
are much reduced, and the resulting filaments are smooth. 

While Hahn et al. (2013) present an elegant alterna¬ 
tive to A-body simulations, their method is numerically 
expensive. Since they are required to track the folding of 
the phase sheet, at the centres of dark matter halos where 
there are many such foliations they formally require an ever- 
increasing number of simulation elements (Hahn Sz Angulo, 
2015); without such refinement, unphysical behaviour occurs 
in high density regions. By contrast, a key strength of the 
‘standard’ A-body method is that, since the equations of 
motion are derived from a Hamiltonian, the time-averaged 
expectation value of the energy of a particle will be correct 
even if its orbital phase is wrong^ (e.g. Dehnen Sz Read, 
2011). This likely explains the success of the A-body method 
even at rather modest A. For example, Sellwood (2006) find 
convergence for their disc simulations already with A ~ 10^, 
despite earlier calculations suggesting some ~ 10® particles 
would be required to properly resolve resonances (Weinberg 
& Katz, 2007). 

The above motivates considering whether the classic A- 
body method cannot be improved. Two recent works have 
attempted to ‘repair’ A-body simulations in post-processing 
by pruning spurious halos. Schneider et al. (2013) propose 
fitting a power law to the artificial halos and subtracting 
them away, taking advantage of the fact that spurious ha¬ 
los are more prevalent in over-dense regions. By contrast, 
Lovell et al. (2014) suggest an algorithm where subhalos are 
removed from the mass function if: (i) their ‘protohalos’ are 
highly flattened; and/or (ii) the subhalos are below a mass 
cut; and/or (hi) the subhalos are not present in a higher 
resolution simulation of the same halo. In this paper, we 
consider instead a modified force softening criterion. This is 
designed to improve the anisotropic force errors that are at 
the root of the problem, leading to a more faithful A-body 
method in the first place. 

This paper is organised as follows. In §2, we briefly re¬ 
view different strategies in the literature for force softening 
and we present our new Novel form of Adaptive softening - 
Nov A - designed to minimise spurious two-body effects. In 
§3, we describe our NovA algorithm in detail and its imple¬ 
mentation in the RAMSES code. In §4, we compare NovA 
to standard RAMSES for an asymmetric plane wave test 
and 0.2 keV WDM simulations. We focus here on present¬ 
ing the first results from NovA for the density field; mass 
function; and dark matter halo density profiles. A detailed 
analysis of halo formation below the WDM ‘free-streaming’ 
mass (see §4.2.1); filament stability; and obtaining new con¬ 
straints on the temperature of dark matter using NovA will 
be presented in forthcoming publications. Finally, in §6 we 
present our conclusions. 


^ This is only strictly true if a symplectic time integrator is used. 
The Leapfrog integrator typically employed in cosmological sim¬ 
ulations is symplectic, but only for fixed timesteps (e.g. Dehnen 
& Read, 2011). 
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2 FORCE SOFTENING 

Since spurious halos in WDM simulations appear to result 
from anisotropic force errors, this suggests that a good place 
to begin in improving the Wbody method is to take a crit¬ 
ical look at how the force softening e is chosen^. Ideally, e 
should be as small as possible such that the maximum pos¬ 
sible force resolution is obtained for a given numerical cost. 
However, too small and spurious forces will creep in, po¬ 
tentially spoiling numerical convergence. We consider three 
different force softening strategies here: 

(i) Minimising two-hody effects (Power): A popular rule- 
of-thumb, that has been carefully tested on cold dark matter 
(CDM) simulations, follows from ensuring that two body 
forces are small as compared to the mean field (Power et ah, 
2003a): 


Gm^ 1 GMm r 

- -e Oi — — 

e2 a‘^r‘^ ^ 


( 2 ) 


where a = 4 is an empirically derived parameter (Power 
et ah, 2003a); N — M/m is the number of particles inside 
some characteristic radius r (Power et al. 2003a use the 
virial radius r 2 oo); and m is the dark matter particle mass. 


(ii) Minimising force errors (Dehnen): An alternative ap¬ 
proach is to minimise errors coming from biased forces that 
occur if e is too large, and noise that occurs if e is too small 
(Dehnen, 2001). This leads to a well-defined optimal force 
softening that depends on the particular gravitational po¬ 
tential being simulated (and the choice of softening kernel). 
For small e and large A", Dehnen (2001) derive an analytic 
estimate for Plummer softening that scales as: 


e (x N valid for e r; N ^ 10^ (3) 

(hi) Minimising scatter between an ensemble of N-body 
realisations (Romeo): Finally, following earlier work by 
Melott et al. (1997) and Splinter et al. (1998), Romeo et al. 
(2008) take a different approach. They run ensembles of the 
same cosmological A-body simulation, varying only the ran¬ 
dom number seed and the force softening. They argue that 
the scatter in results (as measured by various metrics like 
the halo mass function) has a term that is physical (cosmic 
variance), and a term that is numerical (discreteness noise). 
The optimum force softening should minimise the discrete¬ 
ness noise and therefore should minimise the scatter between 
different realisations of the same simulation. Using a novel 
wavelet analysis, they empirically derive: 

e - 2d (4) 

where d is the mean inter-particle spacing. 


^ Note that it is equally important to select an appropriate 
timestep for the particles (e.g. Power et al., 2003a; Zemp et al., 
2007). However, this is true both in the standard A-body method 
and in the ‘folding phase sheet’ model of Hahn et al. (2013). This 
suggests that either the choice of e is more crucial than the choice 
of timestep, or that the timestep criteria typically used in A-body 
simulations (e.g. Dehnen & Read, 2011) are adequate. 


While each of the above approaches to force softening is 
conceptually different, they all point to a rather similar con¬ 
clusion: e must be adaptive, varying both in space and time: 
e = e(x, t). To see this, it is instructive to consider a simple 
toy model where the dark matter is spherically-distributed 
with a power-law density profile: 


peer ^ A(< r) oc (5) 

where A(< r) is the cumulative number of particles within 
r. (This equation is only strictly valid near the centre of 
dark matter halos.) For this toy model, our three criteria 
give rather different results, but all point towards e being 
some function of radius r and therefore of the local density: 


e oc p R > 0 


(6) 


where: 




( 7 — 1)727 Power 

0 . 73(7 — 3)77 Dehnen 
173 Romeo 


There are several interesting points to note from equation 6 . 
Firstly, notice that the Power criterion actually amounts 
to fixed softening if 7 = 1 , as is the case for the centres of 
CDM halos (Dubinski Sz Carlberg, 1991). This may explain 
why fixed softening simulations have performed so surpris¬ 
ingly well despite the natural expectation that e ought to be 
adaptive. If we are to adapt e, however, the Power crite¬ 
rion becomes potentially pathological. For 7 < 1, r < 0 and 
in shallow dark matter cusps or cores, the softening would 
actually increase with density. The Dehnen criterion fares 
better in this regard, being well behaved for all 7 < 3; how¬ 
ever, more work is required to generalise it to larger radii 
where the softening will be large and equation 3 is then no 
longer valid. For these reasons, we consider only the Romeo 
criterion from here on. 

The Romeo softening has (perhaps inadvertently) been 
extensively explored in the literature. Mesh-based methods 
like the RAMSES code (Teyssier, 2002 ) tie the softening to 
the local cell size which is naturally adaptive. Cells are split 
if they have greater than Ac particles, typically chosen to be 
Ac ~ 8 to achieve, on average, one particle per cell after cell 
refinement. This amounts to a scheme where e oc l7Piodi’ 
exactly as in the Romeo force softening. Similar schemes 
have also been explored in Tree A-body codes. There, since 
the equations of motion are derived from a Hamiltonian, it 
is possible to craft a density-adaptive e method that is man¬ 
ifestly conservative (Price & Monaghan, 2007). lannuzzi & 
Dolag (2011) have recently implemented this in the Gadget 
code (Springel, 2005), finding that it leads to results in ex¬ 
cellent agreement with the fixed e case, while giving greater 
resolution for similar numerical cost (see also Kawata et al. 
2013). Their results suggest that with or without the con¬ 
servative correction terms, the halo mass function converges; 
however, without the correction there is a substantial sup¬ 
pression of low-mass halos mass as compared to both the 
conservation-corrected and fixed softening simulations. In 
Appendix A, we show that such conservative corrections do 
not solve ‘spurious halo’ problem in WDM simulations. In 
fact, since the correction terms appear as a purely attractive 



4 Alexander Hobbs, Justin I. Read, Oscar Agertz, Francesca lannuzzi, Chris Power 


Ax 



Figure 1. A schematic view of the need for a modified adaptive 
force softening for A-body simulations. The simulation begins 
with the distribution locally very close to isotropic (left), with 
mean interparticle spacing Ax ~ Ay. However, as collapse pro¬ 
ceeds first along the shortest axis (in this case the x axis), we 
quickly move to a situation where Ax Ay (right). Standard 
adaptive softening schemes adapt purely on the local density. In 
this case, we would pick a softening e oc Ax = making the 
softening too small in the y direction. This could lead to spuri¬ 
ous clumping along the filament. In NovA we measure the local 
anisotropy and do not refine e if the anisotropy is too high. In 
the example pictured, this would set our softening proportional 
to the longest local axis of the collapse - in this case ey oc Ay. 
The softening remains isotropic, but is more conservative than 
standard schemes in regions of high anisotropy. 


force that points along the density gradient, they make the 
spurious halo problem worse. 

In this paper, we present a Novel form of Adaptive 
softening - NovA - designed to minimise spurious two-body 
effects. Like the Romeo softening, we tie the softening to the 
local interparticle spacing e ~ 2d. However, for the first time 
we account for the fact that in cosmological simulations, 
collapse is expected to be locally anisotropic (e.g. Zeldovich, 
1978). Since collapse proceeds most rapidly along the short 
axis of the flow, in the early stages of collapse there will 
always be three interparticle spacings aligned along the short 
c, intermediate h, and long a axis (see Figure 1). If we adapt e 
purely on density, this amounts to an optimistic criterion e ~ 
2c that actually violates the Romeo criterion along the long 
axis, leading to potentially large spurious two-body forces. 
Instead, our new NovA method allows e to be adapted on 
density only if the collapse is sufficiently isotropic. We show 
that this simple change to the A-body algorithm prevents 
the formation of ‘spurious’ halos along filaments. (Note that, 
while in this paper we adapt on density, in principle NovA 
can be applied to any adaptive softening scheme that obeys 
equation 6 . We defer such generalisations to future work.) 


gravitational potential is calculated by solving the Poisson 
equation using the multi-grid method (Guillet & Teyssier, 
2011 ) for all refinement levels. 

Note that while we describe RAMSES as a ‘standard 
A-body’ code, it actually dihers from most Tree A-body 
codes in an important respect: the softening e is automat¬ 
ically adapted according to the Romeo criterion if the re¬ 
finement strategy is based on reaching a critical number of 
particles per cell. Since this is true for all adaptive mesh re¬ 
finement schemes in the literature to date, however, we still 
refer to this as ‘standard’. We compare and contrast RAM¬ 
SES with some Gadget Tree A-body simulations that use 
both fixed and adaptive softening in Appendix A. 


3.2 The NovA algorithm 

In this section, we describe our Novel Adaptive force soft¬ 
ening algorithm: NovA. This is a modified cell splitting cri¬ 
terion implemented in the RAMSES code. Normally, cells 
are split if they contain more than some critical number of 
particles in a cell: Aceii > Ac. Here, we add an additional 
criterion that the cell is sufficiently isotropic as measured by 
its moment of inertia tensor: 


iVij + zfj) ruj -Xij Vij nij -Xij Zij ruj 

-Xij Uij rrij [xij + zfj) rrij -yij Zij rrij 

-Xij Zij rrij -Vij Zij rrij (xfj + yfj) rrij 

( 7 ) 

where xij = Xi — Xj, yij = yi — yj, Zij = Zi — Zj. 

We compute the eigenvalues of the matrix Ip. Ai, A 2 , A 3 , 
which are sorted so that Ai > A 2 > A 3 , and take the ratio 
Qi = A 1 /A 3 to be a measure of the (spatial) anisotropy in 
the particle distribution. 

Gells are split if Aceii > Ac and qi < r/, where 77 is a 
parameter that controls the amount of anisotropy allowed 
for splitting to occur. The ehect of this is shown schemati¬ 
cally in Eigure 1. By refining only where the particle distri¬ 
bution is locally isotropic, NovA ehectively picks the most 
conservative local force softening. As a result, the soften¬ 
ing in anisotropic regions is always at least twice the longest 
inter-particle spacing in a cell: e ^ 2 max[d]. This means that 
the softening is somewhat overestimated with respect to the 
short axis. However, too large softening only ahects the com¬ 
putational efficiency (since for the same particle number we 
have reduced force resolution); whereas too small softening 
can - through two body ehects - be much more problematic. 
An alternative way to think of the algorithm is that it does 
not refine unless collapse is occurring along all three axes. 
This typically reduces refinement in filamentary or elongated 
structures. 




i=i 


3 NUMERICS 

3.1 The RAMSES ‘standard A-body’ code 

We carry out cosmological A-body simulations using 
the Adaptive Mesh Refinement (AMR) code RAMSES 
(Teyssier, 2002). The collisionless dark matter dynamics are 
evolved using the particle-mesh technique (see e.g. Hockney 
& Eastwood, 1988), with gravitational accelerations com¬ 
puted from the gravitational potential on the mesh. The 


3.3 The choice of 77 and Ac 

In the limit A ^ oc, we would ideally have an anisotropy 
bound of 77 = 1 - i.e. cell splitting is allowed only for purely 
isotropic cells. However, in practice noise in the particle dis¬ 
tribution makes it undesirable to set 77 = 1 exactly. Here, 
we choose as default Ac = 32 which is chosen to ensure that 
there are always enough particles in a cell that li can be 
reliably measured (Ac = 32 ensures at least ~ 3 particles 
per spatial dimension); and we select 77 = 1.08. The latter 
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Nc 


ri 

8 

0.748626 

1.374 

16 

0.304557 

1.152 

32 

0.169258 

1.084 

64 

0.102809 

1.051 

128 

0.0661424 

1.033 

256 

0.0448105 

1.022 

512 

0.0301317 

1.015 


Table 1. Using cell particle noise to choose the anisotropy pa¬ 
rameter r] for a given number of particles per cell Nq. The columns 
show Nc’t the variance in anisotropy parameter for random reali¬ 
sations of a uniform particle distribution aq] and our choice of 77 
derived from this analysis: 77 = 1 + O.SfJg. Our approach has two 
desirable properties: (i) 77 is set by the noise level for a cell; (ii) 
as a result, 77 naturally shrinks with Nq. The grey row marks our 
default choice of parameters. 

number is chosen by drawing 32 particles at random from 
a uniform density distribution and calculating the distribu¬ 
tion function of q. We choose 77 = 1 + 0.5crg where aq is the 
variance in q for this random sampling; similar results for a 
range of Nc are reported in Table 1. This has two desirable 
properties: (i) 77 is set by the noise level for a cell; (ii) as 
a result, 77 naturally shrinks with Nc. Note that the above 
implies that varying 77 with the number of particles in a cell, 
or with the refinement level may give improved performance 
over the fixed 77 scheme we explore here; such improvements 
are beyond the scope of this present work. We explore the 
effect of varying A^ceii and 77 in § 4 . 3 . 

3.4 Numerical performance 

As with any new numerical algorithm, we would be remiss 
not to discuss its numerical cost. We find that the relative 
performance of NovA to RAMSES depends on the precise 
choice of parameters and problem setup. Our default NovA 
scheme is ^10 times faster than RAMSES with iVc = 8 . 
However, this is simply because NovA refines less. When 
comparing with standard RAMSES with Nc — 32 (that re¬ 
fines only one level deeper than NovA; see Table 2), NovA is 
almost the same speed. If we compare NovA and RAMSES 
for simulations where NovA refines similarly to RAMSES, 
NovA is about ~ 20% slower. 


4 RESULTS 

4.1 The asymmetric plane wave (Valinia) test 

We set up initial conditions for an asymmetric plane wave 
test as in Valinia et al. (1997) and Hahn et al. (2013). This 
simple 2D test allows us to make a check on two-body dis¬ 
creteness effects without running a full cosmological volume. 
The plane wave is setup along the x-direction with the fol¬ 
lowing (sinusoidal) phase perturbation in the y-direction: 

= 0COS 

where kp = 27^ jL^ ka = Itt/L and Ca = 0.2. L refers to 
the size of the simulation box. 0 sets the value of expansion 
factorat which the first shell crossing occurs - this is set to 



Figure 2 . N = 5123 run for the asymmetrically-perturbed plane 
wave (Valinia) test, comparing standard RAMSES (top) and our 
default NovA method (bottom). The left panels show the particle 
distribution; the right the AMR refinement map. Notice that in 
RAMSES, the filaments break up into regularly-spaced clumps. 
This occurs because the standard cell-splitting criterion refines 
on the filaments (see top right panel). By contrast, in NovA no 
such refinement occurs (see bottom right panel) and the filaments 
are smooth. Both algorithms refine on the bound structures that 
form at nodes, capturing the same caustic/shell-like structures 
reported in Hahn et al. (2013). 


ttc = 1/7.7 ^ 0.13. The initial particle positions and veloci¬ 
ties were obtained by applying the Zel’dovich approximation 
(Zel’dovich, 1970) to an unperturbed regular Cartesian lat¬ 
tice. 

The plane wave is allowed to evolve under the pure grav¬ 
itational potential of the particles up until a = 1. The results 
for our default choice of Nc = 32 and 77 = 1.08 are shown 
in Eigure 2. The left panels show the particle distribution; 
the right the AMR refinement map. Notice that in RAM¬ 
SES, the filaments break up into regularly spaced clumps. 
This occurs because the standard cell splitting criterion re¬ 
fines on the filaments (see top right panel). By contrast, in 
NovA no such refinement occurs (see bottom right panel) 
and the filaments are smooth. Both algorithms refine on the 
bound structures that form at nodes, capturing the same 
caustic/shell-like structures reported in Hahn et al. (2013). 
Since the bound lumps move to lower mass and smaller spac¬ 
ing with resolution (they are non-convergent), NovA gives 
a more faithful simulation of the correct physics. NovA - 
unlike the standard RAMSES V-body implementation - 
converges much more rapidly with resolution. We will show 
this more quantitively with warm dark matter (WDM) sim¬ 
ulations, next. 


, 7 

X-\- ea-r77 COS kaV 

kl 


( 8 ) 
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4.2 Warm dark matter simulations 

4 . 2.1 Initial conditions and simulation analysis 

The Warm Dark Matter (WDM) simulations were set up 
as in Power et ah (2003b). We used cosmological parame¬ 
ters Qo = 0.27, Qa = 0.73, h = 0.705 and ag = 0.81 at 
z = 0 (Komatsu et ah, 2011). Initial conditions were cre¬ 
ated by generating a statistical realization of a Gaussian 
random density field in Fourier space, with variance given 
by the linear matter power spectrum, and the Zel’dovich ap¬ 
proximation used to compute initial particle positions and 
velocities. A CDM power spectrum was obtained by con¬ 
volving the primordial power spectrum P(k) oc with 

the transfer function appropriate for our chosen set of cos¬ 
mological parameters, computed using the Boltzmann code 
CAMB (see Lewis et ah, 2000). The power spectrum for the 
WDM model was then obtained a la Bode et al. (2001b), 
by filtering the CDM power spectrum with an additional 
transfer function of the form 

/ pWDM /7 X \ 1/2 

’"'“w = ) =[1+ (“'•■)"]" (9) 

where a is a function of the WDM particle mass (eq. 9 in 
Bode et al. (2001b)), k is the wavenumber and u = 1.2 is a 
numerical constant. 

We chose a WDM thermal relic mass of = 0.2 keV. 
Following Schneider et al. (2012), this corresponds to an 
effective Tree-streaming’ scale of: 


A 


eff 

fs 



/ h \ 

(oTr) 


( 10 ) 

which for = 0.25; h = 0.73; and = 0.2 gives Xff = 
0.308 Mpc//i. And a Tree-streaming mass scale’: 



( 11 ) 


where p is the mean background density of the Universe 
(that is a function of cosmology and redshift z; p(z = 0) = 
277.3/i^ M© kpc“^; Peacock 1999). For purely linear collapse 
with no mode-coupling or fragmentation, no halos should 
form below Mfg. At reshift z = 0, and assuming the above 
cosmological parameters and WDM thermal relic mass, we 
have Mfs = 2.26 x 10^ M©. 

A second scale of interest is the length scale at which the 
WDM transfer function is reduced by half: the ‘half-mode’ 
length: 


Ahm 13 . 93 Af (12) 

which also has an associated mass scale, the ‘half-mode 
mass’: 


Mhm 2.7 X lO^Mfs (13) 

This is the mass scale at which we expect the WDM mass 
function to noticeably deviate from the CDM case. 

It is not clear if halos should exist below Mpg in WDM 
structure formation. Angulo et al. (2013) use their new ‘fold¬ 
ing phase sheet’ method to argue that there are no halos be¬ 


low Mfs (other than substructure halos that originate from 
halos more massive than Mfg). However, this result relies on 
some manual pruning of the halo mass function - required 
due to errors in the halo finding algorithm. We verify that 
halo finding in WDM is indeed a thorny issue, and discuss 
Nov A results for halos below Mfg in §4.2.2. 

Note that we deliberately choose a small = 0.2 keV 
even though such a low thermal relic mass is already ruled 
out by constraints from the Lyman-a forest (e.g. Viel et ah, 
2013). The reason for this is that it ensures that WDM ef¬ 
fects will appear at large mass, making it computationally 
efficient to test our methodology (c.f. Hahn et ah, 2013). 
We will present Nov A simulations of particle masses close 
to the current observational constraints (and comparisons 
with data) in future work. A full list of all simulations run 
and their parameters is given in Table 2. 


4 . 2.2 The density field 

Figure 3 shows a full box view of the projected dark matter 
density field in RAMSES and NovA. Bound halos are iden¬ 
tified using AHF (see §4.2.1); these are overplotted as red 
filled circles, with a size proportional to their virial radii. In 
the RAMSES simulation (left panels), the filaments break 
up into many small halos. These ‘spurious’ halos have al¬ 
ready been shown to be numerical artefacts (e.g., Wang & 
White, 2007) and it is encouraging that they are gone in 
NovA (bottom panels). Once these structures are removed, 
NovA does a good job of capturing the caustics, fine shells, 
and criss-crossing filaments that surround galaxies. 

Where caustic and shell structures overlap, AHE often 
identifies ‘halos’, yet it is not clear if such structures are 
real or simply transient. We call these ‘caustic’ halos. Their 
existence - if real - is important as it implies that halos can 
indeed form below the free-streaming mass Mfg. A detailed 
exploration of this requires improving the halo finder and 
studying these structures carefully as a function of time. 
We will explore this in detail in a forthcoming publication. 


4 . 2.3 The halo mass function and convergence 

In Eigure 4, we compare the AHE cumulative halo mass 
function for RAMSES (blue) with NovA (red) at three 
numerical resolutions: N — 256^ (dashed), 512^ (dotted) 
and 1024^ (solid) particles. Overplotted is the cumulative 
mass function for the equivalent cold dark matter simula¬ 
tion (black), and the free streaming mass Mfg at redshift 
^ = 0 (dotted line). 

The RAMSES simulations show very poor numerical 
convergence, with a prominent upturn that shifts extremely 
slowly to lower mass as the resolution is increased. It is in¬ 
teresting that our results for mass function convergence in 
RAMSES are somewhat worse than reported in Wang & 
White (2007). This likely owes to the fact that we use a 
different halo finder; we find that switching off the ‘unbind’ 
feature in AHE leads to very different behaviour, illustrating 
how sensitive results for WDM simulations are to the choice 
of halo finder and its chosen settings. 

By contrast, NovA shows much better numerical con¬ 
vergence. The mass function rises slowly at the low-mass end 
with resolution, while remaining unchanged at high mass. 
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Figure 3. A comparison of RAMSES and NovA for a 0.2 keV WDM simulation with N = 512^ particles. The top panels show the 
full 50Mpc/h box; the bottom panels highlight a zoomed in region, as marked. Halos identified using AHF are marked by the red filled 
circles; their size is proportional to their virial radii. Notice that in RAMSES, many small and regularly spaced halos - ‘spurious halos’ 
- form along the filaments; in NovA these are gone. Notice further that in the zoom panel for NovA, there are several cases of halos 
identified by AHF that may not correspond to genuine bound structures that could host galaxies. We highlight two of these ‘caustic’ 
halos as examples. Some of these halos likely owe to overlapping caustics in the WDM density field and would not host galaxies; others, 
however, may be genuine non-linear structures that form at an overlap between caustics or as filaments intersect. This latter possibility 
is very interesting as it would imply that halos can form below the free streaming mass Mfg in WDM. We will explore this further in a 
forthcoming paper. 


Below the free streaming scale, we find nearly an order of 
magnitude fewer halos in NovA than in RAMSES. We do 
however find a tendency for the lowest resolutions to over¬ 
suppress halos at intermediate mass, although this goes away 
quickly with increasing resolution. 


Our chosen halo finder AHF has been extensively tested 
on CDM simulations (Knebe et ah, 2011), but has not been 
used on ‘spurious halo free’ WDM simulations before. As 
discussed in §4.2.2 (and see Figure 3), for the simulations 
we present here, AHF almost certainly misidentifies some 
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Table 2. The cosmological simulations and their parameters. The columns show from left to right: the simulation label; the particle 
resolution; the minimum number of particles per cell the anisotropy parameter 77 (see §3); the dark matter particle mass rupart; and 
the maximum refinement level reached by that simulation. 


Label 

Resolution 

Nc 

77 

TUpart (IO^^Mq/Zi) 

Max. refinement level 

RAMSES-256 

256=^ 

32 

none 

7.1 X 10-^ 

14 

RAMSES-512 

5123 

32 

none 

8.9 X 10-3 

16 

RAMSES-1024 

10243 

32 

none 

1.1 X 10-3 

17 

NovA-256 

2563 

32 

1.08 

7.1 X 10-^ 

13 

NovA-512 

512® 

32 

1.08 

8.9 X 10-3 

15 

NovA-1024 

10243 

32 

1.08 

1.1 X 10-3 

16 

NovA-256-NcI28 

256® 

128 

1.03 

7.1 X 10-2 

7 

NovA-5I2-NcI28 

5123 

128 

1.03 

8.9 X 10-3 

8 

NovA-I024-NcI28 

10243 

128 

1.03 

1.1 X 10-3 

14 


features as halos. For this reason, we defer a more careful 
analysis of halos below the free streaming mass Mfg to future 
work. We discuss how our results compare with other recent 
determinations in the literature in §5. 

Note that our NovA algorithm does lead to increas¬ 
ing refinement with resolution. This is shown in the final 
column of Table 2, where we list the maximum refinement 
level reached for all WDM simulations in this paper. The 
NovA simulations typically reach ~ 1 refinement level less 
than standard RAMSES, but nonetheless they do continue 
to refine with increasing resolution. We discuss the effect 
of varying 77 on the maximum refinement reached and on 
numerical convergence in §4.3 and §5.1. 


4 . 2.4 Halo density profiles 

In Figure 5, we show dark matter halo density profiles in 
RAMSES and NovA for the 256^ and 512^ simulations for 
both low mass 10^^ M©) and high mass (~ 10^^ M©) 
halos. It has already been reported previously in the lit¬ 
erature that simply forbidding refinement will reduce spuri¬ 
ous two-body effects in collisionless iV-body simulations (e.g. 
Melott et ah, 1997; Splinter et ah, 1998; Hahn et ah, 2013). 
However, reducing refinement everywhere also leads to the 
centres of halos - where galaxies actually reside - becom¬ 
ing unresolved. In NovA, we attempt to obtain the best of 
both worlds as the algorithm leads to derefinement in highly- 
anisotropic regions while still refining on halo centres (see 
Figure 2, bottom right panel). For this reason, the NovA 
halos typically have a density profile in excellent agreement 
with the RAMSES simulation over all of the resolutions 
studied here. The lowest mass mass halos in NovA are shal¬ 
lower than their RAMSES counterparts reflecting the lower 
refinement level reached. The effect diminishes with resolu¬ 
tion, however, demonstrating that NovA is convergent. 


4.3 Varying Nc and 77 : Should some filaments 
fragment after all? 

In this section, we study the effect of varying the minimum 
number of particles in a cell Nc and the anisotropy param¬ 
eter 77 . For this paper, we use throughout 77 = 1 + O.barj as 
outlined in §3. We defer a detailed analysis of 77 - in par¬ 
ticular allowing a time or spatially varying 77 - to future 
work. 



log [M (M3,,)] 

Figure 4. Dark matter cumulative halo mass functions in RAM¬ 
SES (blue) and NovA (red) for a 0.2 keV WDM simulation. In 
both cases, three different resolutions are marked: N = 256^ 
(dashed); N = 512^ (dotted); and N = 1024^ (solid) particles. 
Overlaid is the curve expected for the equivalent cold dark matter 
simulation (black). Notice also that the standard RAMSES sim¬ 
ulations converge very slowly with increasing resolution, showing 
a characteristic upturn at low mass that only strengthens as the 
resolution increases. By contrast, NovA converges rapidly ‘from 
below’. At A = 1024^ particles, the number of halos at low mass 
in NovA is suppressed with respect to standard RAMSES by 
nearly an order of magnitude. Finally, notice that even in NovA, 
the cumulative mass function does not reach a plateau below the 
free-streaming mass scale: there are significant numbers of ha¬ 
los below Mfs (vertical dotted line; see §4.2.1). See the text for 
further discussion of this. 

In Figure 6 , we show a zoom-in on a particularly mas¬ 
sive filament that can be seen in the top middle of the full 
simulation box shown in Figure 3. We focus on this partic¬ 
ular filament because we find that it is very hard to avoid it 
fragmenting. At low resolution (256^; left panels), the fila¬ 
ment is completely smooth. However, for our default choice 
of Nc = 32, already at A = 512^, the filament begins to 
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Figure 5. Dark matter density profiles in a 0.2 keV WDM simulation for example halos with virial masses lO^'^Msun (upper lines) and 
~ lO^^Msun (lower lines) in RAMSES (blue) and NovA (red) in the 256^ (left) and 512^ (right) resolutions. The profiles correspond to 
the same halos at each resolution in each method. 


fragment (top middle panel), as does another massive fila¬ 
ment that connects to the large halo at the edge of the box 
(see yellow circles). As we raise the numerical resolution N 
at fixed Nc and 77 , the filament fragments further (top right 
panel). If we raise Nc and lower 77 according to our ‘noise 
criteria’ (§3; bottom panels) then at iV = 512^ the filament 
becomes once again smooth (bottom middle panel). How¬ 
ever, even for these NovA parameters, raising the resolution 
to N = 1024^ results in the filament fragmenting (bottom 
right panel). Indeed, the filament breaks up into structure in 
any version of NovA where we permit refinement within the 
filament. Interestingly, the largest structures that form both 
in this filament and the one to its top left appear always in 
the same place regardless of our choice of NovA parame¬ 
ters. These structures are marked by the yellow circles (yel¬ 
low is chosen to avoid confusion with AHF halos in previous 
plots that are marked in red). The fact that such structures 
are challenging to avoid and yet also appear always at the 
same locations suggests that they may be physically cor¬ 
rect. Such behaviour is certainly rather different from the 
regularly spaced fragments that form along filaments in the 
standard RAMSES simulations (Figure 3). We discuss this 
further in §5. 


5 DISCUSSION 

5.1 Numerical Convergence in NovA 

It is challenging to quantitatively test convergence in NovA 
because of uncertainties in the halo finding. For example, 
Wang & White (2007) using Friends-of-Friends (FoF) report 
a definite, if slow, shift of the upturn in the halo mass func¬ 
tion to lower mass with increasing resolution. Using AHF, 
we find also a shift in the upturn but it is substantially 
slower (see Figure 4). Taken at face value, it is not at all 
clear that RAMSES will converge on the correct solution 
with increasing A" - at least not when using AHF. 

By contrast, convergence in NovA seems healthier. At 


low mass, the mass function is suppressed causing a ‘con¬ 
vergence from below’ with increasing N. At low resolution, 
there is a clear over-suppression in halos even at high mass 
(at N = 256^ RAMSES outperforms NovA). But with in¬ 
creasing A, the situation rapidly improves. We find a gently 
rising mass function with few low mass halos, but nonethe¬ 
less some. Again, owing to difhculties with halo finding, we 
cannot yet determine whether these low mass halos are real; 
a fault of the halo finder; or evidence that NovA requires 
further improvement. We defer a careful analysis of this to 
a future work where we will improve on the halo finder and 
study the time evolution of halos. 


5.2 Comparison with recent work in the literature 

Despite the nearly an order of magnitude suppression in 
spurious halos in NovA, at first sight the mass function ap¬ 
pears to rise substantially more steeply towards low mass 
than that reported recently in Angulo et al. (2013); using 
their new ‘folding phase sheet’ methodology, they find no 
halos below Mfg. However, there are three confounding fac¬ 
tors that make a comparison difhcult. Firstly, Angulo et al. 
(2013) use a ‘Friends of Friends’ (FoF) halo finder similarly 
to Wang & White (2007), whereas we use AHF (that in¬ 
cludes, for example an ‘unbinding’ procedure that discards 
halos that are not gravitationally self-bound). Secondly, they 
throw out all halos with overlapping virial radii. This has the 
effect that all substructure halos are removed, providing a 
better comparison with semi-analytic theories for structure 
formation like Press-Schechter (Press & Schechter, 1974). 
However, such a method will also remove half of all bi¬ 
nary systems. Finally, they further prune their mass function 
through a visual inspection of halos. This is necessary due 
to the FoF algorithm picking up many false positives. We 
find fewer false positives when using AHF, but many of the 
issues they report with halo finding in WDM resonate with 
our findings here. As pointed out in §4.2.2 (and see Figure 
3), we do see many suspicious structures identified by AHF 
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Figure 6 . Zoom-in on a particularly massive filament (taken from the top middle of the full simulation box; see Figure 3) for varying 
Nc] T]] and numerical resolution N, as marked on each panel. We focus on this particular filament because we find that it is very hard to 
avoid it fragmenting. At low resolution (256^; left panels), the filament is completely smooth. However, for our default choice of Nc = 32, 
already at A” = 512^, the filament begins to fragment (top middle panel), as does another massive filament that connects to the large 
halo at the edge of the box (see yellow circles). As we raise the numerical resolution N at fixed Nc and 77 , the filament fragments further 
(top right panel). If we raise Nc and lower 77 according to our ‘noise criteria’ (§3; bottom panels) then at A = 512^ the filament becomes 
once again smooth (bottom middle panel). However, even for these NovA parameters, raising the resolution to A = 1024^ results in the 
filament fragmenting (bottom right panel). Interestingly, the largest structures that form in filaments appear in the same place regardless 
of our choice of NovA parameters; these are marked by the yellow circles. 


as bound halos that are unlikely to host galaxies. Given the 
difficulty of halo finding in WDM, we defer a more careful 
analysis of halos below Mfg to a future work where we will 
look critically at the time evolution of halos in NovA and 
attempt to improve on existing halo finding algorithms that 
have all been tuned to work well only on CDM simulations 
(e.g. Knebe et ah, 2011). 

5.3 Filaments that physically fragment? 

In crafting a scheme that avoids spurious fragmentation of 
filaments, we should be mindful of not throwing the baby 
out with the proverbial bathwater. Are we really sure that 
no filament should fragment, or that all such fragments re¬ 
ported in A-body simulations are spurious? Our results in 
§4.3 strongly suggest that at least some filaments may be un¬ 
stable to physical fragmentation. From a theoretical point 
of view, this should perhaps not be surprising. It has long 
been known that infinite self-gravitating cylinders are un¬ 
stable to fragmentation (Ostriker, 1964; Fridman & Poli- 
achenko, 1984). However, cosmological filaments are both 


finite and expanding^ and given sufficient expansion, self- 
gravitating cylinders become unconditionally stable (Schnei¬ 
der & Moore, 2011). We will explore this further in fu¬ 
ture work, but the issue is an important one: if filaments 
are physically unstable then halos can collapse below Mfg 
and isolated small galaxies can be expected to exist even in 
warm/hot dark matter cosmologies. 


6 CONCLUSIONS 

We have introduced a Novel form of Adaptive softening 
(NovA) for collisionless A-body simulations, implemented 
in the RAMSES adaptive mesh refinement code. In RAM¬ 
SES - that we refer to as a ‘standard A-body method’ - 
cells are only split if they contain more than eight particles 
(a mass refinement criterion). We introduced an additional 
criterion that the cell be sufficiently isotropic, as measured 
by the ratio of the maximum to minimum eigenvalues of its 
moment of inertia tensor: 77 = Amax/Amin. In this way, col¬ 
lapse is only refined if it occurs along all three axes, ensuring 
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that the softening e is always of order twice the largest inter¬ 
particle spacing in a cell. This more conservative force soft¬ 
ening criterion was designed to minimise spurious two-body 
effects, while maintaining high force resolution in collapsed 
regions of the flow. 

We tested NovA using an antisymmetric perturbed 
plane wave collapse (‘Valinia’ test) before applying it to 
warm dark matter (WDM) simulations. Our key results are 
as follows: 

• We used the Valinia test to show that - unlike the stan¬ 
dard Wbody method (RAMSES) - NovA produces no nu¬ 
merical fragmentation while still being able to correctly cap¬ 
ture high density features like the fine caustics and shells 
around the collapsing regions. 

• For the WDM simulations, we found that NovA con¬ 
verges significantly more rapidly than RAMSES, producing 
little or no spurious halos on small scales. NovA produces 
nearly an order of magnitude fewer dark matter halos at low 
mass as compared to RAMSES, while still being able to cor¬ 
rectly resolve high density regions at the centres of massive 
halos. 

• Despite the large reduction in low mass halos in NovA, 
we still found halos below the free streaming mass scale Mfg. 
Furthermore, these halos increase in number (albeit slowly) 
as we increase the numerical resolution. Some of these likely 
owe to our halo finder (AHF) incorrectly labelling caustics 
and criss-crossing filaments as halos. Others form as larger 
halos that form above Mfg are tidally stripped. However, 
some isolated low-mass structures appear to be real. Due 
to the difficultly of accurate halo identification in WDM, 
we defer a quantitative analysis of halos below Mfg to a 
forthcoming publication. 

• We highlighted two particularly massive filaments that 
fragment in any version of NovA where refinement is al¬ 
lowed. Since the most massive fragments appear always at 
the same locations, we argue that these may be physical. We 
noted that infinite self-gravitating cylinders are unstable to 
collapse and so particularly massive cosmological filaments 
may be physically unstable. We will explore this further in 
future work, but the issue is an important one: if filaments 
are physically unstable then halos can collapse below Mfg 
and isolated small galaxies can be expected to exist even in 
warm/hot dark matter cosmologies. 

We will use NovA in forthcoming papers to study the issue 
of halo formation below Mfg; filament stability; and to obtain 
new constraints on the temperature of dark matter. 
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APPENDIX A: WHY CONSERVATIVE 
ADAPTIVE SOFTENING DOES NOT SOLVE 
THE PROBLEM OF SPURIOUS HALOS 


One of our earlier ideas for solving the spurious halo problem 
was to use conservative adaptive force softening as originally 
suggested by Price & Monaghan (2007) and implemented re¬ 
cently in the Gadget code (Springel, 2005) by lannuzzi & 
Dolag (2011). In this Appendix, for completeness, we report 
the results of this experiment and explain why it failed to 
produce the desired effect. In fact, such conservative cor¬ 
rections to adaptive force softening make the problem of 
spurious halos even worse! 

As detailed in Price & Monaghan (2007), if the force 
softening varies in space and time then we can still con¬ 
struct a fully conservative V-body method by deriving the 
equations of motion from a discretised Hamiltonian. This 
results in an additional corrective force to the usual V-body 
equation of motion: 


dvi,c 

dt 

where: 


and: 


U. dWijhj 0 dWijhA 

2 ^ ^ dvi dvj 

j L J J 


(Al) 
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dhi dWij(hi) 
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(A3) 


where G is Newton’s gravitational constant; ruj is the mass 
of particle j; Wij is a positive definite spherical smoothing 
kernel; hj is the smoothing length of particle j, here equated 
with the softening hj = ej; rj is the position of particle j; 
and (t)ij = (j)(\ri — rj\) is related to the gravitational potential 
between particle pairs (Price & Monaghan, 2007). 
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Figure Al. Gumulative mass function for Gadget adaptive soft¬ 
ening runs, with (magenta) and without (brown) the conservative 
correction term, compared to the fixed softening reference run 
(black). Notice that the upturn is more pronounced when the 
correction term is used. These runs used N = 128^ particles. 


The key thing to note from equation Al is that, since 
the kernel is positive definite, and similar terms 

are negative definite; the correction terms Qi are always of 
order unity; while the ^i terms are also negative definite. This 
means that the force correction in equation Al is negative 
definite and will lead always to an increased gravitational 
force. This force will point to leading order along the density 
gradient: 


^ipi = '"^rnjViWijihi) (A4) 

j 

This is the trouble with the conservative correction terms. 
When the flow becomes highly anisotropic, the correction 
terms will act to increase the force along the short axis lead¬ 
ing to even more artificial clumping (see Figure 1). 

In Figure Al, we show the cumulative mass function for 
a 50Mpc/h 0.2 keV WDM simulation run using the standard 
Gadget iV-body code (black); with adaptive force soften¬ 
ing (brown); and with conservative adaptive force softening 
(magenta). For this comparison, we used a kernel neighbour 
number of iVneigh = 60 since this was shown to be optimal 
in lannuzzi & Dolag (2011). Notice that the sharp upturn in 
the mass function (present in all simulations) is more promi¬ 
nent in the simulation with the conservative correction term 
(magenta) than even in the fixed softening case (black). 















